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We study the clustering of passive, non-interacting particles moving under the influence of a 
fluctuating field and random noise, in one dimension. The fluctuating field in our case is provided 
by a surface governed by the Kardar-Parisi-Zhang (KPZ) equation and the sliding particles follow 
the local surface slope. As the KPZ equation can be mapped to the noisy Burgers equation, the 
problem translates to that of passive scalars in a Burgers fluid. We study the case of particles 
moving in the same direction as the surface, equivalent to advection in fluid language. Monte- 
Carlo simulations on a discrete lattice model reveal extreme clustering of the passive particles. The 
resulting Strong Clustering State is defined using the scaling properties of the two point density- 
density correlation function. Our simulations show that the state is robust against changing the 
ratio of update speeds of the surface and particles. In the equilibrium limit of a stationary surface 
and finite noise, one obtains the Sinai model for random walkers on a random landscape. In this 
limit, we obtain analytic results which allow closed form expressions to be found for the quantities 
of interest. Surprisingly, these results for the equilibrium problem show good agreement with the 
results in the non-equilibrium regime. 
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I. INTRODUCTION 



The coupling of two or more driven diffusive systems 
can give rise to intricate and interesting behavior, 
and this class of problems has attracted much recent 
attention. Models of diverse phenomena, such as growth 
of binary films Jl|, motion of stuck and flowing grains 
in a sandpile |2|, sedimentation of colloidal crystals 
and the flow of passive scalars like ink or dye in 
fluids 0, involve two interacting fields. In this 
paper, we concentrate on semiautonomously coupled 
systems — these are systems in which one field evolves 
independently and drives the second field. Apart from 
being driven by the independent field, the passive field 
is also subject to noise, and the combination of driving 
and diffusion gives rise to interesting behavior. Our 
aim in this paper is to understand and characterize the 
steady state of a passive field of this kind. 

The passive scalar problem is of considerable interest 
in the area of fluid mechanics and has been well studied, 
see 0, for reviews. Apart from numerical studies, 
considerable understanding has been gained by analyzing 
the Kraichnan model where the velocity field of a 
fluid is replaced by a correlated Gaussian velocity field. 
Typical examples of passive scalars such as dye particles 
or a temperature field advected by a stirred fluid bring 
to mind pictures of spreading and mixing caused by the 



combined effect of fluid advection and diffusion. On the 
other hand, if the fluid is compressible, or if the inertia 
of the scalars cannot be neglected 0, the scalars may 
cluster rather than spread out. It has been argued that 
there is a phase transition as a function of the com- 
pressibility of the fluid — at large compressibilities, the 
particle trajectories implode, while they explode in the 
incompressible or slightly compressible case |g- It is the 
highly compressible case which is of interest in this paper. 

Specifically, we study and characterize the steady 
state properties of passive, non-interacting particles 
sli ding on a fluctuating surface and subject to noise 
9, 101. The surface is the autonomously evolving field 
and the particles slide downwards along the local slope. 
We consider a surface evolving according to the Kardar- 
Parisi-Zhang (KPZ) equation. This equation can be 
mapped to the well known Burgers equation with noise, 
which describes a compressible fluid. Thus the problem 
of sliding passive particles on a fluctuating surface maps 
to the problem of passive scalars in a compressible fluid. 
We are interested in characterizing the steady state of 
this problem, first posed and studied by Drossel and 
Kardar in ||. Using Monte-Carlo simulations of a solid 
on solid model and analyzing the number of particles in 
a given bin as a function of bin size, they showed that 
there is clustering of particles. However their analysis 
does not involve the scaling with system size, which 
as we will see below, is one of the most important 
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characteristics of the system. We find that the two point 
density-density correlation function is a scaling function 
of r and L (r is the separation and L is the system size) 
and that the scaling function diverges at small r/L. 
The divergence indicates formation of clusters while the 
scaling of r with L implies that the clusters are typically 
separated from each other by a distance that scales with 
the system size. A brief account of some of our our 
results has appeared in [Icj . 

Scaling of the density-density correlation function 
with system size has also been observed in the related 
problem of particles with a hard core interaction, 
sliding under gravity on a KPZ surface 0, 0, H^j |. 
However, the correlation function in this case has a cusp 
singularity as r/L — > 0, in contrast to the divergence 
that we find for noninteracting particles. Thus, while 
clustering and strong fluctuations are seen in both, the 
nature of the steady states is different in the two cases. 
In our case, clustering causes a vanishing fraction of sites 
to be occupied in the noninteracting case, whereas hard 
core interactions force the occupancy of a finite fraction. 
In the latter case, there are analogies to customary 
phase ordered states, with the important difference that 
there are strong fluctuations in the thermodynamic 
limit, leading to the appellation Fluctuation Dominated 
Phase Ordering (FDPO) states. The terminology Strong 
Clustering States is reserved for the sorts of nonequilib- 
rium states that are found with noninteracting particles 
- a key feature being the divergent scaling function 
describing the two point correlation function. 

In the problem defined above, there are two time 
scales involved, one associated with the surface evolution 
and the other with particle motion. We define lo as the 
ratio of the surface to the particle update rates. While 
we see interesting variations in the characteristics of the 
system under change of this parameter, the particular 
limit of lo — ► is of special importance. There is a slight 
subtlety here as the limit lo — > does not commute 
with the thermodynamic limit L — > oo. If we consider 
taking lo — ► first and then approach large system size 
(L — * oo), we obtain a state in which the surface is 
stationary and the particles move on it under the effect 
of noise. In this limit of a stationary surface, we obtain 
an equilibrium problem. This is the well known known 
Sinai model which describes random walkers in a random 
medium. We will discuss this limit further below. Now 
consider taking the large system size limit first and 
then approach lo — 0; this describes a system in which 
particles move extremely fast compared to the evolution 
of the local landscape. This leads to the particles 
settling quickly into local valleys, and staying there till a 
new valley evolves. We thus see a non-equilibrium SCS 
state here, but with the features that the probability 
of finding a large cluster of particles on a single site 



is strongly enhanced. We call this limiting state the 
Extreme-strong clustering state (ESCS) (Fig. 0. The 
opposite limit shown in Fig. ^ is the lo — > oo limit 
where the surface moves much faster than the particles. 
Because of this very fast movement, the particles do not 
get time to "feel" the valleys and they behave as nearly 
free random walkers. 

The equilibrium limit (lo — > followed by L — > oo) 
coincides with the Sinai model describing random 
walkers in a random medium [Tij ]. This problem can be 
analyzed analytically by mapping it to a supersymmetric 
quantum mechanics problem and we are able to 

obtain closed form answers for the two quantities of 
interest — the two point correlation function G(r, L) 
and the probability distribution function of finding n 
particles on a site P(n,L). Surprisingly, we find that 
not only do these results show similar scaling behavior 
as the numerical results for lo = 1 (nonequilibrium 
regime) but also the analytic scaling function describes 
the numerical data very well. The only free parameter 
in this equilibrium problem is the temperature and we 
choose it to fit our numerical data for the nonequilibrium 
system. Interestingly, the effective temperature seems to 
depend on the quantity under study. 

The KPZ equation contains a quadratic term which 
breaks the up-down symmetry, thus one can have 
different behavior of the passive scalars depending 
on whether the surface is moving downwards (in the 
direction of the particles, corresponding to advection 
in fluid language) or upwards (against the particles, 
or anti-advection in fluid language). In this paper, 
we will consider only the case of advection. One can 
also consider dropping the nonlinear term itself; this 
leads to the Edwards- Wilkinson (EW) equation for 
surface growth. The problems of KPZ anti-advection 
and passive sliders on an Edwards Wilkinson surface 
are interesting in themselves and will be addressed in a 
subsequent paper [l6| . 

Apart form the static quantities studied above, one 
can also study the d yna mic properties of the system. 
Bohr and Pikovsky |l7l | and Chin have studied 

a similar model with the difference that they do not 
consider noise acting on particles. In the absence of 
noise, all the particles coalesce and ultimately form a 
single cluster in steady state, very different from the 
strongly fluctuating, distributed particle state under 
study here. References 01 an d 01 study the process 
of coalescence in time. Further, they find that the RMS 
displacement for a given particle increases in time t as 
i 1 / 2 , where z is equal to the dynamic exponent of the 
surface, indicating that the particles have a tendency 
to follow the valleys of the surface. Drossel and Kardar 
[9j have studied the RMS displacement in the same 
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its motion is governed by 



dx r 
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FIG. 1: Change in state with change in u. In the u) —* 
limit one gets different kinds of states depending on how one 
approaches it. The uj — » oo is the free particle limit. The 
arc shoes that there is a similarity between the results of the 
equilibrium Sinai limit and the non-equilibrium SCS at u) = 1. 




FIG. 2: Schematic diagram of the surface and non-interacting 
particles sliding on top of it. Arrows show possible surface and 
particles moves. 



problem in the presence of noise and observe the same 
behavior. We confirm this result in our simulations and 
observe that the variation of u> does not change the result. 

The arrangement of this paper is as follows. In Sec- 
tion II, we will describe the problem in terms of contin- 
uum equations and then describe a discrete lattice model 
which mimics these equations at large length and time 
scales. We have used this model to study the problem 
via Monte Carlo simulations. Section III describes re- 
sults of our numerical simulations. We start with results 
on the various static quantities in the steady state and 
define the SCS. We also report on the dynamic quantities 
and the effect on steady state properties of varying the 
parameter w. Section IV describes our analytic results 
for the equilibrium Sinai limit of a static surface and the 
surprising connection with results for the noncquilibrium 
problem of KPZ/Burgers advection. 



II. DESCRIPTION OF THE PROBLEM 

The evolution of the one-dimensional interface is de- 
scribed by the KPZ equation [J^ 



dh 
~di 



d 2 h 
dx 2 



A dh. 



2 + Ch(x,t). 



(1) 



Here h is the height field and C,h is a Gaussian white noise 
satisfying ({ h (x, tX h (x', t')) = 2D h 6 d {x-x')d(t~t'). For 
the passive scalars, if the m th particle is at position x m , 



dt 



dh 

dx 



Cm(t) 



(2) 



where the white noise Cm(t) represents the randomiz- 
ing effect of temperature, and satisfies (Cro(*)Cm(*')) = 
2nS(t — t'). Equation J5J is an overdamped Langevin 
equation of a particle in a potential h(x, t) that is also 
fluctuating, with a determining the speed of sliding. In 
the limit when h(x,t) = h(x) is static, a set of noninter- 
acting particles, at late times would reach the equilibrium 
Boltzmann state with particle density ~ e -0 ,l ( x ) _ On the 
other hand, when h(x,t) is time dependent, the system 
eventually settles into a strongly nonequilibrium steady 
state. The transformation v = —dh/dx maps Eq. (JTJ 
(with A = 1) to the Burgers equation which describes a 
compressible fluid with local velocity v 



dv 
dt 



d 2 v , d( h (x,t) 



dx 2 



dx 



(3) 



The above equation describes a compressible fluid be- 
cause it does not have the pressure term, which is present 
in the Navier Stokes equation. The transformed Eq. 
describes passive scalar particles advected by the Burgers 
fluid 



dt 



aV \x^ +Cm(t) 



(4) 



The ratio a/A > corresponds to advection (particles 
moving with the flow), the case of interest in this paper, 
while a/A < corresponds to anti-advection (particles 
moving against the flow). 

Rather than analyzing the coupled Eqs. (JIJ and 
or equivalently Eqs. © and Q directly, we study a 
lattice model which is expected to have similar behavior 
at large length and time scales. The model consists of a 
flexible, one-dimensional lattice in which particles reside 
on sites, while the links or bonds between successive 
lattice sites are also dynamical variables which denote 
local slopes of the surface. The total number of sites 
is L. Each link takes either of the values +1 (upward 
slope — ► /) or —1 (downward slope — * \). The rules for 
surface evolution are : choose a site at random, and if it 
is on a local hill (— > /\), change the local hill to a local 
valley(^ \/) (Fig. |2J). After every N s surface moves, we 
perform N p particle updates according to the following 
rule : we choose a particle at random and move it one 
step downward with probability (1 + K)/2 or upward 
with probability (1 — K)/2. The parameter K ranges 
from 1 (particles totally following the surface slope) to 
(particles moving independently of the surface). In 
our simulations, we update the surface and particles at 
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independent sites, reflecting the independence of the 
noises Ch(x, t) and ( m (t) [2(J]. The ratio u = N s /N p 
controls the relative time scales of the surface evolution 
and particle movement. In particular, the limit ui — ► 
corresponds to the adiabatic limit of the problem where 
particles move on a static surface and the steady state 
is the thermal equilibrium state. 

To see how the lattice model described above describes 
a KPZ surface, consider the mapping of the above model 
to the well known asymmetric simple exclusion process 
(ASEP): consider an up slope to be a particle on a lattice 
and a down slope to be an empty space (hole). The flip- 
ping of a hill to a valley then corresponds to the motion 
of a particle (exchange of particle and hole). A coarse 
grained description of the ASEP leads to the KPZ equa- 
tion [2l]]. The continuum description of the ASEP, ob- 
tained by coarse graining over regions which are large 
enough to contain many sites, involves the density of 
particles p(x) and the local current J(x). These are con- 
nected through the continuity equation 



dp 

at 



aj 

d.r 



The local current can be written as 



t/ \ dp • < \ 

J(x) = -V— + j(p) + T) 



(5) 



(6) 



where v is the particle diffusion constant, 77 is a Gaussian 
noise variable and j(p) is the systematic contribution to 
the current associated with the local density p. Using 
the expression for the bulk ASEP with density p for j, 
we have 



j(p) = (P - l)p( l - P) 



(T) 



where p and q are the particle hopping probabilities to 
the right and left respectively, with our one-step model 
corresponding to p = 1 and q = 0. 



the signs of the constant term and A are opposite. Thus 
a downward moving surface (corresponding to p > q) 
has positive A. The constant term can be eliminated by 
the boost h — > h — i(p — q)t, but its sign is important 
in determining the overall direction of motion of the 
surface. The case (a/ A) > which is of interest to us 
thus corresponds to the lattice model in which particles 
move in the same direction as the overall surface motion. 

The parameters to and K defined in the lattice model 
are connected to the continuum equations as follows. In 
the limit of a stationary surface, we achieve equilibrium 
and the particles settle into in a Boltzmann state with 
particle density ~ e -P h ( x ) ^ here h(x) is the surface height 
profile and [3 is the inverse temperature. (3 is related to 
K by /3 = In ^y^r^ and to the parameters a and k in 
Eq. © by (3 = a/n. Thus 



K 



_ ^ 

■Mj k _|_ y 



(10) 



The parameter u> cannot be written simply in terms 
of the parameters in the continuum equations, because 
it modifies Eq. as we now show, w is the ratio of 
the update speeds or equivalently the time between suc- 
cessive updates of the particles (At p ) and surface (At s ). 
The noises C,h(x,t) and ( m (t) in Eqs. and © can 

be written as yj ~j^£h( x -> t) and y^f-CmW respectively. 

Here Ch(x, t) is noise of O(l), uncorrelated in time, white 
in space while Cm(t) is uncorrelated noise of 0(1). The 
factors of \J~-^[ in the terms indicate that the strength of 
the noise depends on how frequently noise impulses are 
given to the particles; the square root arises from the ran- 
dom nature of these impulses. Thus the change in height 
(Ah) in time At s and the distance traveled (Ax m ) in 
time Atp are respectively - 



Ah = At s [, 



8 2 h 
dx 2 



A dh 2 

2 W J 



VAt s D h ( h (x,t) (11) 



Since we identify the presence (absence) of a particle 
in the lattice model with an up (down) slope, we may 
write 



P : 



I., dh, 
2 (1 + ^ 

Using Eqs. ©,0 and © in Eq. © leads to 



dh 1 d 2 h 1, s,dh. 2 

m = -2 {p - q) + "d^ + 2 {p - q){ d^ ) 



(8) 



(9) 



which is the KPZ equation (Eq. ©) with an additional 
constant term, and A = (p — q) and — —rj. Note that 



Ax n 



At r , 



dh 



dx n 



+ y/At pK ((t) (12) 



We now identify At s and At p with the Monte-Carlo time 
step St as At s = N S S and At p = N P S. We can thus re- 
place At s by uj.Atp and take it to be the natural contin- 
uous time. We thus get 



dh 
dt 



, d 2 h A,on. 9 , _ 

- = + o (— ) ] + v^aOM) (13) 



A ,9/iv 2 i 
2 ( &? 



dh 



dt 



dx r 



<m(t) 



(14) 
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We can see that the u> dependence in the above equation 
cannot be removed by a simple rescaling of the param- 
eters of the equation. Eq. Q is recovered as a special 
case of Eq. I|13fl on setting u> = 1. 



III. NUMERICAL RESULTS 
Two Point Density Density Correlation Function 

We start with the simplest case u> = K = 1; surface 
updates are attempted as frequently as particle updates, 
and both particles and surface always move only down- 
wards. In our simulations, we work with N = L, where 
N is the total number of particles and there are L sites 
in the lattice. The two point density-density correlation 
function is defined as G(r, L) — (njnj+ r ) l, where rii is 
the number of particles at site i. Fig. [3] shows the scaling 
collapse of numerical data for various system sizes (L) 
which strongly suggests that for r > 0, the scaling form 



(15) 



is valid with 6 ~ 1/2. The scaling function Y(y) has a 
power law divergence Y(y) ~ y~ u as y — > 0, with v close 
to 3/2. The data for r = points to G(0, L) ~ L. 

This numerical result matches with an exact result 
of Dcrrida et. al. jl^l for a slightly different model. 
As we have seen in the previous section, the single step 
model which we use for Monte-Carlo simulations, can be 
mapped on to an asymmetric simple exclusion process 
(ASEP). The particles/holes in the ASEP map to the 
up/down slopes in our model and the flipping of a hill 
to a valley is equivalent to swapping a particle with a 
hole. In (2^, apart from particles and holes, a third 
species called the second-class particles are introduced 
which act as holes for the particles and particles for 
the holes. When translated to the surface language, 
these second class particles behave like the sliders in our 
model, with the difference that they are not passive: 
there is no surface evolution at a site where second-class 
particles reside. The effect of non-passivity is relatively 
unimportant for KPZ advection-like dynamics of the 
surface, as particles mostly reside on stable local valleys 
while surface evolution occurs at local hilltops. More- 
over, if the number of second class particles is small, 
the probability of the rare event where they affect the 
dynamics of local hills goes down even further. With 
only two such particles in the full lattice, probability 
p(r) that they are at a distance r, is proportional to the 
two point correlation function G{r,L). The exact result 
|22| p(r) r~j ^jry matches very well with our prediction 
for the same quantity, p(r) = -^G(r, L) <~ 7372. 




FIG. 3: The inset shows G(r, L) versus r for different values 
of L. The main plot shows the scaling collapse when r is 
scaled with L and G(r,L) with 1/L ' 5 . The dashed, straight 
line shows y ~ x~ 1,5 . The lattice sizes for both plots are L= 
256 (*), 512 (x), 1024 (+). 



The result for G(r, L) also allows us to calculate 
the quantity N(l,L) first defined in 9J; the lattice is 
divided into L/l bins of size I and we ask for the number 
N(l,L) of particles in the same bin as a randomly 
chosen particle. N(l,L) is a good measure of clustering 
- if N(l, L) rises linearly with I, one concludes that the 
particles are fairly spread out, while if N(l,L) saturates 
or shows a decreasing slope, one concludes that particles 
are clustered. N(l, L) is related to the two point 
correlation function through N(l, L) = G(r, L)dr, 
using which we obtain N(l, L) ~ c\L{\ — C2l~ v+1 ). This 
form fits the numerical result for iV better (Fig. |3J than 
the ^-independent form of 0. 



Probability Density of Occupancy 

Another quantity of primary interest is the probabil- 
ity P(n, L) that a given site is occupied by n particles. 
For n > 0, this quantity shows a scaling with the total 
number of particles, which in turn is proportional to the 
system size L. We have (see Fig. 



P(n,L) 



— f(- 



(16) 



with 5 = 1. The scaling function f(y) seems to fit 
well to a power law y~ 7 with 7 ~ 1.15 (Fig. EJ, 
though as we shall see in Section IV, the small y 
behavior may follow y^ 1 lny. We can use the scal- 
ing form in the above equation to calculate G(0, L), 
(n 2 ) = G(0,L) = J L n 2 P(n,L)dn - L s = L, which, as 
we have seen above, is borne out independently by the 
numerics. Numerical data for P(0, L) (which is not a 
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FIG. 4: The inset shows Nil, L) versus bin size I for different 
system sizes (L). The main plot shows N(l,L) scaled with 
L versus bin size I. The curve shows C\L{1 — C2l~ u+1 ) with 
ci = 1 and C2 = 0.72. The straight line shows N(l,L) — L, 
the form predicted in [IJ. The lattice sizes for both plots are 
L= 256 (*), 512 (x), 1024 (+). 
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FIG. 5: The inset shows P(n, L) versus n for different val- 
ues of L. The main plot shows L 2 P(n,L) versus n/L. The 
straight line shows y ~ a; -1 ' 15 . The lattice sizes are L= 256 
(*), 512 (x), 1024 (+). 
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FIG. 6: The inset shows V(P(l,L)) versus P(l,L) for differ- 
ent values of L, L= 256 (*), 512 (x), 1024 (+). The main plot 
shows V(P(1, L)) versus P(l, L) for different number of aver- 
aging times, t= T/10 (+), T (*), T * 5 (x) where T=30,000. 

L — ► oo. Let us ask for the probability density function 
V(P(n, L)) describing the values taken on by P(n, L). As 
seen in (Fig.^J, this distribution does not change when 
we increase the averaging time (main figure) or the length 
(inset). Thus r P(P(n,L)) approaches a distribution with 
a finite width in the thermodynamic limit rather than a 
delta function. This clearly indicates that there are large 
fluctuations in the system which do not damp out in the 
thermodynamic limit. Large fluctuations, which do not 
decrease with increasing system size, are also a feature of 
the FDPO state for particles with a hard core interaction 

[HEG3. 

Results on Dynamics 

The root mean square (RMS) displacement R(t) = 
((x(t) — z (O)) 2 ) 1 / 2 of a tagged particle has been studied 
earlier pHl"8j|. R(t) is found to obey the scaling form 



part of the scaling function in Eq. (|16fl ) shows that the 
number of occupied sites N occ = (1 — P(0, L))L varies as 
with (f) ~ 0.23, though the effective exponent seems 
to decrease systematically with increasing system size L. 



Fluctuations 

To evaluate the fluctuations of the quantity P(n,L) 
about its mean, we evaluated the standard deviation — 
AP(n, L) = ^(P(n,L) 2 ) - (P(n,L)) 2 where the brack- 
ets denote an average over time. We find that this quan- 
tity does not decrease even in the thermodynamic limit 



R{t) = L^hU-\ (17) 

where h(y) ~ y 1 ^, with z = 3/2 for small y. The 
requirement that R(t) has to be independent of L in 
the limit L — > oo leads to x — 1- The value of z 
above is the same as the dynamic exponent of the KPZ 
surface. The dynamic exponent z s of a surface carries 
information about the time scale of evolution of valleys 
and hills; the landscape evolves under surface evolution 
and valleys/hills of breadth L' are typically replaced by 
hills/ valleys in time of order L' Zs . Thus the observation 
z = z s suggests that the particles follow the valley 
movement. 
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We have also evaluated the autocorrelation function 
G(t,L) = (iii(0)ni(t)) l and find that it scales with the 
system size as 

G(t,L)^Y(±y (18) 

Again, z = z s = 3/2, reaffirming our conclusion that par- 
ticles tend to follow valleys. The scaling function shows 
a power law behavior Y(y) ~ with -0~2/3asy^O. 



Relations Between the Exponents 

The exponents defined in the above sections can be 
connected to each other by simple relations using scaling 
analysis. For instance, S, v and are related by 

5 = v- 9 (19) 

This can be proved by substituting the scaling form of 
Eq. H3) and G(0,L) = J^n 2 P{n,L)dn ~ L 5 in the 

equation L G(?*, L)dr — L; the last equation can be ob- 
tained by using the definition of G(r, L) and using N = L. 
We can also relate (j>, S and 7 by 

= 6{j - 2) + 1 (20) 

which can be derived using the normalization condition 
P(n, L)dr = 1 and then substituting for P(0, L) 
and the scaling form of Eq. I|16|) . Our results from 
simulations are consistent with these relations. 

The following picture of the steady state emerges from 
our results. The scaling of the probability distribution 
P(n,L) as n/L and the vanishing of the probability of 
finding an occupied site (= N occ /L) suggest that a large 
number of particles (often of the order of system size) 
aggregate on a few sites. The scaling of the two-point 
density-density correlation function with L implies 
that the particles are distributed over distances of the 
order of L, while the divergence of the scaling function 
indicates clustering of large-mass aggregates. Thus the 
evidence points to a state where the particles form a 
few, dense clusters composed of a small number of large 
mas aggregates and these clusters are separated on the 
scale of system size. We choose to call this state as the 
Strong Clustering State (SCS). The divergence at origin 
of the two-point density-density correlation function as 
function of the separation scaled by the system size, 
is its hallmark. The information we get from results 
on dynamics is that the particles have a tendency to 
follow the surface. This is brought out by the fact that 
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FIG. 7: Scaled probability distribution P(n, L) for ui — 
1/2, 1,2 (if = 1). The line is a fit to Eq. g3 with = 2.3. 
The lattice sizes are L= 512 (o, x, □), 1024 (■, +, *). 
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FIG. 8: The main plot shows the scaled two point correlation 
function for ui = 2, (K = 1), we see deviation from scaling 
at small r/L. The inset shows a plot of G(r,L)/L versus r. 
The straight line shows depicts the power law y ~ a; -1 ' 5 . The 
lattice sizes are L= 256 (+),512 (x), 1024 (*). 

the scaling exponent describing the RMS displacement 
comes out to be equal to the dynamic exponent of the 
KPZ surface. 



Variation of u and K 

To see how the system behaves when we change the rel- 
ative speeds of the surface and particle evolution, we vary 
the parameter uj = N s /N p (N s and N p being respectively 
the number of successive surface and particle update at- 
tempts) in the range 1/4 < to < 4. When u> < 1 (parti- 
cles faster than the surface), we regain the scaling form 
of Eq. I|15|) for the two point correlation function. The 
scaling function also diverges with the same exponent. 
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While the probability distribution for occupancy P(n, L) 
shows similar scaling with system size as Eq. <|16ll . the 
scaling function f(y) shows a new feature — it develops 
a peak at large n (Fig. 0) . This peak at large n indicates 
that the probability of finding nearly all the particles at 
a single site is substantial. A heuristic argument for the 
appearance of this peak is now given. Consider a config- 
uration in which a large number of particles (nearly equal 
to the total number of particles) reside in a local valley. 
When this valley is replaced by another one nearby un- 
der surface dynamics , all the particles tend to move to 
the new one. If the number of particle updates is greater 
than surface updates, there is a substantial probability 
that all the particles are able to move to the new valley 
before it is itself replaced by another one. Thus there 
is a significant probability of the initial cluster surviving 
intact for a fair amount of time. Numerically, we also 
find that 




FIG. 9: Two point scaled density correlation G(r, L) function 
(advection) for u) = 1, K = 0.75 (top curve), 1 (bottom). 
The line is a plot of Eq. 173H with j3 = 4. The lattice sizes are 
L = 1024 (□, x), 512 (*, +). 



P(n = N) 
P(n = N - 1) 



1 

u> 



(21) 



For u> > 1, the particles settle down slowly in val- 
leys and T sur f T part where T sur f and T part are respec- 
tively the times between successive surface and particle 
updates. Though T sur f ^> T part ; for large enough L, the 
survival time of the largest valley ~ r sur fL z is always 
greater than the particle sliding time ~ T part L. Thus we 
expect that particles will lag behind the freshly formed 
valleys of small sizes but would manage to cluster in the 
larger and deeper valleys, which survive long enough. We 
thus expect a clustering of particles and scaling to hold 
beyond a crossover length scale (r c (uj)). We can estimate 
the crossover length by equating the time scales of surface 
and particle rearrangements — T sur fr z (uj) ~ T part r c (uj), 
which yields r c (ui) ~ oj^ 1 . Using z = 3/2, we have 
r c ~ u> 2 . Numerical simulation results are shown in 
Fig. 03 which shows that the data deviates from scaling 
and power law behavior at small r, due to the crossover 
effect. The data suggests that 



G(r,L) ~~ jjY 



(z>< 



(22) 



As we can see from Fig. [S] (main graph), the curve 
flattens out at small values of r, so for y < 1 (r < r c (u>)), 
the function g(Y) in the equation above should follow 
g(y) ~ y 15 while it should go to a constant for y > 1. 
We can determine r c (u>) from G(r,L) by separating out 
the r dependent part; if we scale G(r, L) by L, we obtain 
the quantity -przg(^rj^))- We can now determine r c (uS) 
as the value or r where the scaled data starts deviating 
from the power law behavior r" 15 . From Fig. EI (inset) 
r c (oj = 2) ~ 10. A similar exercise for oj = 3 leads to 
r c (uj = 3) ~ 20. A clean determination of r c {ui) for 



u> > 3 requires data for very large values of system size, 
beyond the scope of our computational capabilities. 

The probability distribution P(n, L) continues to 
show the same scaling form (Eq. (|16(l ) for u > 1, but 
the scaling function /(y) in this case dips at large values 
of y (Fig. EJ in contrast to the peak seen for oj < 1. 
The exponent z describing the RMS displacement of 
particles remains unchanged under a change in u>, again 
indicating that particles follow the movement of valleys 
on the large scale. 

The other parameter of interest is K, defined in Sec- 
tion II — when we make a particle update, we move the 
particle downhill with probability (1 + K)/2 and uphill 
with probability (1 — K)/2. So far we have discussed the 
results for the case K = 1, where particles always move 
downhill. Decrease in K reduces the average downhill 
speed of particles, while the valley evolution rates are 
unaffected. Thus decreasing K causes an effect similar 
to increasing w and a crossover length scale is introduced. 
The particles lag behind the freshly formed local valleys 
but settle down in the deeper, longer surviving global 
valleys. The numerical results again guide us to the form 



(23) 



for the correlation function. Analogous to u> > 1 case, 
we have extracted r c from the numerical data. We find 
r c (K — 0.75) ~ 10. Values of K lower than 0.75 require 
data for system sizes that are beyond our computational 
limitations. 
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IV. ADIABATIC, EQUILIBRIUM LIMIT 



The Exact Distribution of the Probability Density 



To approach the problem analytically we take the ex- 
treme limit of uj — > 0. In this limit, the surface is sta- 
tionary. The particles move on this static surface under 
the effect of noise and the problem becomes a well known 
equilibrium problem - the Sinai model for random 
walkers on a random landscape. It is well known that 
for the KPZ surface in one-dimension, the distribution of 
heights h(r) in the stationary state is described by 



Prob[{/i(r)}] ex exp 



2D h 



dhlr') 
dp 



dr' 



(24) 



Thus, any stationary configuration can be thought of as 
the trace of a random walker in space evolving via the 
equation, dh(r)/dr = £(r) where the white noise £(r) has 
zero mean and is delta correlated, (£( r )£( r ')) = S(r — r'). 
We impose periodic boundary conditions as for the lat- 
tice model, without loss of generality - h(0) — h(L) = 0. 

The passive particles moving on top of this surface as, 
we remember, move according to Eq. Since this is 
an equilibrium situation, (( m (t)£ m (t')} — 2n5(t — t') = 
2KsT5(t — t') where T is the temperature and Kb is 
the Boltzmann constant. Since the particles are non- 
interacting, we can deal with a single particle and instead 
of the number of particles n r at a site r, we consider the 
probability p(r)dr that the particle will be located be- 
tween r and r+dr. In the long time time, the particle will 
reach its thermal equilibrium in the potential h(r) and 
will be distributed according to the Gibbs-Boltzmann dis- 
tribution, 



p(r) 



-/3h(r) 



(25) 



where Z = J L dre~^ h ^ is the partition function. Note 
that p(r) in Eq. (|25|l depends on the realization of the 
potential {h(r)} and varies from one realization to an- 
other. Our goal would be to compute the distribution 
of p(r) over different realization of the random potential 
h(r) drawn from the distribution in Eq. I)24[l. Note that 
the distribution of h(r) in Eq. (|24Jl is invariant under 
the transformation h(r) — > —h{r). In other words, the 
equilibrium density p(r) defined in Eq. (|25() will have 
the same distribution if one changes the sign of h(r) in 
Eq. H25fl . For later simplicity, we will make this trans- 
formation now and replace p(r) instead by the following 
definition 



p(r) 



a0Hr) 



(26) 



where the transformed partition function is now given by 
Z = /'drew. 



Our strategy would be first to compute the n-th mo- 
ment of the random variable p(r) in Eq. (|26() . It follows 
from Eq. $3$ - 



P n {r) 



nf3h(r) 



r(n) 



dyy n -h 



-yZ+n0h(r) 



(27) 



where we have used the identity dyy n 1 e yZ = 
T(n)/Z n to rewrite the factor l/Z n . Here T(n) is the 
standard Gamma function. Next, we make a further 
change of variable in Eq. H27|) by writing y — f3 2 e l3u /2. 
Note that as y varies from to oo, the new dummy vari- 
able u varies from — oo to oo. Making this substitution 
in Eq. |23> we get, 

p n (r) = b n J duexp[-^|^ dx e W*)+«)j 

+nj3(h(r) + u)} (28) 

where we have used the explicit expression of the par- 
tition function, Z = dre 13 ' 1 ^. The constant b n = 
P 2n+1 /[2 n T(n)]. We are now ready to average the ex- 
pression in Eq. I|28() over the disorder, i.e., all possible 
realizations of the random potential h{x) drawn from the 
distribution in Eq. 124(1 . Taking the average in Eq. 1281) 
(we denote it by an overbar), we get using Eq. 124(1 . 



ph(L)=o 
du / T>h{x) exp[ 

Jh(0)=0 
2 



p n (r) = Ab r , 




1 f dh(x)\ | (3 2 c 0( h ( x)+u) 



2\ dx J 2 

n/3(h(r) + u)] 



(29) 



where the normalization constant A of the path integral 
in Eq. 1(29(1 will be chosen so as to satisfy the sum rule, 

Jo p(r)dr = 1. Next we shift the potential by a constant 
amount u, i.e., we define a new function V(x) — h(x) + u 
for all x that reduces Eq. 1(29(1 to the following expression, 



p n (r) 




I : path integral can be viewed as a quantum mechan- 
ical problem in the following sense. All paths (with the 
measure shown above) starts from V(0) = u and ends 
at V(L) = u. At the fixed point r (where we are trying 
to calculate the density distribution), these paths take a 
value V(r) = V which can vary from — oo to oo. This 
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can be written, using the quantum mechanical bra-ket 
notation, 



oo poo 



p n {r) = Ab n du / dV < u\e- Hr \V > e n0V 

J —oo J — oo 
-H(L-r 



< V\e 



u > 



(31) 



The first bra-ket inside the integral in Eq. (|31[) denotes 
the propagation of paths from the initial value u to V at 
the intermediate point r and the second bra-ket denotes 
the subsequent propagation of the paths from V at r to 
the final value u at L. The Hamiltonian H corresponds 
to the operator H = \ (^) 2 + ^-e 0V{ - x \ Interpreting 
V(x) to be the "position" of a fictitious particle at the 
fictitious "time" x, this operator has a a standard kinetic 
energy term and a potential energy which is exponential 
in the "position" V. The right hand side of Eq. I|31[) can 
be rearranged and simplified as in the following - 



pn( r ) = Ab n I dVe nf}v / du<V\e- H( - L - r 1\u> 

) 

(32) 



<u\e- Hr \V> 



Thus, 



p-n(r) = Ab n dVe npv < V\e- HL \V > (33) 

J — oo 

where we have used the completeness condition, 
/_ du\u >< u\ — I with I being the identity oper- 
ator. At this point, it may be helpful and less confusing 
notationally if we denote the "position" V of the fictitious 
quantum particle by a more friendly notation V = X, 
which will help us thinking more clearly. Thus, the Eq. 
(|33|l then reduces to, 



p n (r)=Ab r , 



dXe npX < X \ £ -HL, X > 



(34) 



To evaluate the matrix element in Eq. (|34|) . we need to 
know the eigenstates and the eigenvalues of the Hamil- 
tonian operator H . It is best to work in the "position" 
basis X. In this basis, the eigenfunctions iPe(X) of H 
satisfies the standard Shrodinger equation, 



2 dX 2 2 6 



tI>e(X) = E^ E (X), (35) 



valid in the range — oo < X < oo. It turns out that this 
Shrodinger equation has no bound state (E < 0) and 
only has scattering states with E > 0. We label these 
positive energy eigenstates by E — /3 2 k 2 /8, where k is 
a continuous label varying from to oo. A negative k 
eigenfunction is the same as the positive k eigenfunction, 
and hence it is counted only once. With this labeling, 
it turns out that the differential equation can be solved 
and one finds that the eigenfunction ip^{X) is given by, 



MX)=a k K lk (2eW 2 ) , 



(36) 



where K u (y) is the modified Bessel function with index v. 
Note that, out of two possible solutions of the differential 
equation, we have chosen the one which does not diverge 
as X — > oo, one of the physical boundary conditions. The 
important question is: how to determine the constant 
afe in Eq. Note that, unlike a bound state, the 

wavefunction ^fc(X) is not normalizablc. To determine 
the constant a k , we examine the asymptotic behavior of 
the wavefunction in the regime X — » — oo. Using the 
asymptotic properties of the Bessel function (when its 
argument we find that 



MX) 



T(ik) 



-ik0X/2 



a ikl3X/2 



I 2sin(i/br)r(l + ik) 

(37) 

On the other hand, in the limit X — > — oo, the 
Schrodinger equation (|35|l reduces to a free problem, 



l d 2 MX) i?k 2 
which allows plane wave solutions of the form, 



MX) 



^ikf3X/2 



"(fc)< 



-ikpX/2 



(38) 



(39) 



where e lfc ^ x / 2 represents the incoming wave from X — 
— oo and e ~ lk x / 2 represents the reflected wave going 
back towards X — — oo with r(k) being the reflection 



coefficient. The amplitude \ f- is chosen such that the 



plane waves ipk(X) = y -j^e tk P x l 2 are properly orthonor- 

malized in the sense that < ^kW k >= S(k — k') where 
S(z) is the Dirac delta function. Comparing Eqs. IpTfjl 
and (|39J) in the regime X —> — oo, we determine the con- 
stant a,k (up to a phase factor), 



/3 

sin(i/c7r)r(l 



ik). 



(40) 



The square of the amplitude |afc| 2 (which is independent 
of the unknown phase factor) is then given by 



\ak\ 



[3k sinh(7rfc) 



(41) 



where we have used the identity, T(l + ik)T((\ — ik) = 
7rfc/sinh(7rfc). Therefore, the eigenstates of the operator 
H are given by \k >, such that H\k >— ^-£-\k > and in 
the X basis, the wavefunction ijjk(X) =< k\X > is given 
(up to a phase factor) by the exact expression 



(42) 



We now go back to Eq. I|34|l where we are ready to eval- 
uate the matrix element < X\e~ HL \X > given the full 
knowledge of the eigenstates of H . Expanding all the 
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kets and bras in the eigenbasis \k > of H, we can rewrite 
Eq. (P3J as follows, 

/OO /"OO 
dX dk < X\k >< k\X > 
-oo JO 

e n(3X e -/3 2 k 2 L/8 

dX\MX)\ 2 e n ^ x . (43) 



= Ab n 



dk i 



The X integral on the right hand side of Eq. (|43[) can be 

expressed in a compact notation, 



dX\MX)\ 2 e nl3x =< k\e nf)jt \k > . (44) 



Substituting the exact form of ipk{X) from Eq. I|42|) . wc 
get 

fc sinh(7rfc) 



< k\e^ x \k >= n J^_l> j n dyy 2n ^K ik {y)K^ k {y). 

(45) 

Fortunately, the integral on the right hand side of Eq. 
(|45J) can be done in closed form [23J and we obtain, 



< k\ 



3 n/3X 



k > = 



fcsinh(Trfc) T 2 {n) 
tt^ 2 "- 1 r(2n) 



r(n — ik)T{n + ik). 



Substituting this matrix element back in Eq. 143|) . 
arrive at our final expression, 



(46) 



p n {r) = A 



dk k sinh(7rfc) 



o 

2 1.2 , 



4tt 2 2" T(2n) 
\T(n-ik)\ 2 e~^ L /\ 
To determine the constant A, we first put n 



(47) 



1 in Eq. 



(|T7|l . Note that p(r) = 1/L by virtue of the probability 
sum rule, J Q L p(r)dr — 1. Taking the disorder average and 
using the translational invariance, one gets p(r) = 1/L. 
Using the identity, T(l + ik)T((l — ik) — irk/s'mh(irk) 
and performing the integral on the right hand side of Eq. 
(|47|l and then demanding that the right hand side must 
equal 1 /L for n = 1, we get 



A = V2nL. 



(48) 



One can also check easily that n — > 0, the right hand side 
of Eq. I|47[) approaches to 1 as it should. In verifying this 
limit, we need to use the fact that T(x) ~ 1/x as x — > 
and also the identity, T(ik)Y(—ik) — 7r/fcsinh(7r£;). Now, 
for n > (strictly), one can make a further simplification 
of the right hand side of Eq. I)47|) . We use the property 
of the Gamma function, T(x + 1) = xT(x), repeatedly 
to write T(n — ik) = (n — 1 — ik)T(ri — 1 — ik) = (n — 
1 - ik)(n - 2 - ik) . . . (1 - ik)T(l - ik). Note that this 
formula, so far, is valid only for integer n > 1. This gives, 
for integer n > 1 

r(n - ik)T(n + ik) = [(n - l) 2 + k 2 }[(n - 2) 2 + k 2 }... 

irk . 



where we have used the identity, T(l + ik)T((l — ik) — 
7rfc/sinh(7rfc). Substituting this expression in Eq. I|47|l wc 
get, for n > 1, 



p n (r) = V2ttL 



p2n+l r(n) 



4tt2™ T(2n 
[(n-2) 2 + fc 2 ]... 



dkk 2 [{n- l) 2 + fc 2 ] 
[1 + k 2 ]e-^ k2L '\ (50) 



Making the change of variable f3 2 k 2 L/8 = z in the in- 
tegral, we finally obtain the following expression for all 
integer n > 1, 



p n (r) 



1 /? 2 



I» 



iv^F 2™- 2 r(2n) 



rfz ( 



,1/2 



2 2 H 




/3 2 L_ 





(n - l) 2 + 



8z 



P 2 L 
(51) 



For example, consider the case n = 2. In this case, the 
formula in Eq. (|51(l gives 



P 2 

p 2 (r) = 



12 
~]PL 



(52) 



which is valid for all L and not just for large L. Note that 
the second term on the right hand side gives a contribu- 
tion which is exactly 1/L 2 . This means that the variance, 

2 

p 2 (r) — p{r) = f3 2 /[12L] for all L. For arbitrary integer 
n > 1, taking the large L limit in Eq. (|51|l we get, as 
L —>■ oo, 



P n (r) -> i 



2n-2 



r 3 (n) 



ln-2 



r(2n) 



(53) 



. Note that even though this expression was derived 
assuming integer n > 1, after obtaining this formula, 
one can analytically continue it for all noninteger n > 0. 
Now, let us denote Prob(p, L) = P(p,L). Then p n (r) = 
J °° p n P(p, L)dp. Note again that the range of p is from 
to oo , since it is a probability density, and not a proba- 
bility. The factor 1 /L on the right hand side of Eq. (|53|) 
suggests that P(p, L) has the following behavior for large 



P(p,L) = I f(p), 
where the function f(y) satisfies the equation, 



y n f(y)dy 



T(2n) 



(54) 



(55) 



To determine f(y) from this equation, we first use the 
identity, T(2n) = 2 2n - 1 r(n)r(?i + 1/2)/^, known as 
the doubling formula for the Gamma function. Next we 
use El, 
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Identifying the right hand side of Eq. (|56|l with the right 
hand side of Eq. upon choosing a — 2/[3 2 , we get 

the exact expression of f(y), 



2 

Wy 



More cleanly, we can then write that for large L, 



■2p 



where the scaling function f'(y) is universal (independent 
of the system parameter 0) and is given by 



f'(y) = —K (y). 
V 



(59) 



This function has the following asymptotic behaviors, 



f'(y) 



i[-ln(j//2) + 0.5772...], 




0, 
oo. 



(60) 



The scaling form in Eq. I|58(l is valid only when m(r) ~ 
L. If m(r) is a number of order O(l) (not as large as L), 
then the scaling breaks down. This fact suggests that the 
correct behavior of the distribution P(p, L) for large L 
actually has two parts, 



P(P,L) 



1 - 



ln 2 (L) 



(3 2 L 



2p 

d 2 



('-!)■ 

(61) 



where f'(y) is given by Eq. I|59|l . This form in Eq. 1|61|) 

is consistent with all the observed facts. For example, if 
one integrates the right hand side, the first term gives 1 — 
ln p2f? (with the convention S(y)dy =1). The second 

term, when integrated, gives ^gff (where we have used 
the small y behavior of f'(y) from Eq. (|60|l and kept only 
the leading order term for large L) which exactly cancels 
the identical factor in the first term to give a total sum 1, 
as it should. On the other hand, for any finite moment of 
order n, the first term does not contribute and only the 
second term contributes to give the result in Eq. J53J). 



The Density-Density Correlation Function 

We now consider the density-density correlation func- 
tion between two points r\ and r 2 at equilibrium. The 
calculation proceeds more or less along the same lines as 
in the previous section. The density-density correlation 
function is defined as 



C(n,r 2 ) = p(ri)p(r 2 ), 



(62) 



which evidently depends only on r = \r% — r 2 \ due to the 
translational invariance. The density p(r) is again given 
by Eq. (gSJl. It follows from Eq. ® that 



(57) p(n)p(r 2 ) 



(58) 



o P[h(n)+h(n)] 



Z 2 



dyye 



-yZ+p[h{ri)+h{r 2 )] 

3 

(63) 

where the partition function, Z = j Q dre^ u ^ and we 
have used the identity 1/Z 2 = f^°dyye~~ Zy . As in 
Section-II, we now make a change of variable in Eq. I|63|l 
by writing y = /3 2 e /3u /2. Then Eq. 163|l becomes, 



p(ri)p(r 2 ) 



/3 5 [°° 

— du exp[ 

^ J —OO 

f3(h{n)+u + h{r 2 )+u)], 




(64) 



where we have used the explicit expression of the par- 
tition function, Z = dre^ r \ Averaging over the 
disorder, we get 



p(n)p(r 2 ) 



B 



4 

1 / dh(x 



h(L)=0 

du I T>h(x) exp[ 

h(0)=0 




1 fdh{x) 
2 V dx J ^2 V dx 

+(3(h{ ri ) + h(r 2 )+2u)} 



+ 2 



(65) 



where the normalization constant B will be determined 
from the condition, J J C{ri,r 2 )dridr2 = 1 (which 

follows from the fact that J Q L p(r)dr — 1). Alternatively, 
one can put r = r 2 — r\ = in the expression for the 
correlation function and then it should be same as p 2 (r) 
already computed in the previous section. 



As before, we next shift the potential, i.e., we define 
V(x) = U (x) + u for all x. The Eq. I|65|l then simplifies, 



p{ri)p(r 2 ) 



B 



du 



V{L)=u 



VV(x) 



V(0) = 









exp 


H 


I dx _ 









§_ p l3V(x) 



1 fdV(x) 
2 1 dx 



+ /3(V(n) + V(r 2 )) 



(66) 



Thus we have again reduced the problem to a path 
integral problem. However, there is a difference in the 
subsequent calculations. This is because, unlike the pre- 
vious calculation, we now have to divide the paths into 3 
parts: (i) paths starting at V(0) = u and propagating up 
to the point r\ where ^(^i) = V\ (note that V\ can vary 
from — oo to oo), (ii) paths starting at r\ with V(r\) = V± 
and propagating up to r 2 with V(r 2 ) — V 2 and (iii) paths 
starting at r 2 with V(r 2 ) = V 2 and propagating up to L 
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where V(L) = u. We have assumed r 2 > r% for conve- 
nience. Using the bra-ket notation, we can then re-write 
Eq. (PI as 



p(n)p{r 2 ) 



D 



ff 



dm 



dVi 



dVo 



To extract the asymptotic behavior for large L, we 
rescale k% y/L — r = x\ and fc 2 VL = X2 in Eq. i|70|) . then 
expand the sinh's and the cosh's for small arguments, 
perform the resulting double integral (which becomes 
simple after the expansion) and finally get for L — > oo 
and r^O, 



m 



' 2 < V 2 \e~ H( - L ^\u > 



(67) 



The Hamiltonian H 
in the previous section. Using J 
(|6T|) can be simplified, 

,/3 5 



I (iD + ^ V{x) is the same as 



du\u X u\ = I, Eq. 



p(ri)p(r 2 ) = B'— I dV x 



where r = r% — T\. Note that Eq. 1|68[) clearly shows 
that C(n,r2, L) = C(r = r 2 — r\, L), as it should due to 
the translational invariance. Furthermore, Eq. Ij68(l also 
shows that that function C (r, L) is symmetric around 
r = L/2, i.e., C{r,L) = C(L - r,L). This last fact is 
expected due to the periodic boundary condition. As 
before, we change to a more friendly notation: V\ = X\ 
and V2 = X 2 , where Xi and X 2 denote the 'positions' of 
the fictitious quantum particle at 'times' r\ and r 2 . With 
this notation, Eq. i|tj8fl reads, 



P{r\)p{ri) 



D 



dX x 



dX 2 



< X 2 \e-^ L -^\X X > 

< X!\e- Ar \X 2 > e ^ Xl+X2 l (69) 



This can be solved to obtain the correlation function 



C{r,L) = B 



256 



dkidk2kik2(kf — k 2 ) 



sinh(7rfci) sinh(7rfc2 



[cosh(7rfci 
exp 



cosh(7rfc2)] 5 



klr) 



(70) 



For r = 0, it is possible to perform the double integral in 
Eq. l|7UI) and one finds that it reduces to the expression 
of p 2 (r) in Eq. Q52J1 of the previous section, provided the 
normalization constant B = \ 2ttL. Thus, the two-point 
density-density correlator is given exactly by Eq. H7U|) 
(with B = \/2itL) and note that this expression is 
valid for all L. This exact expression of the correlation 
function was first derived by Comtet and Texier 0] 
in the context of a localization problem in disordered 
supersymmetric quantum mechanics. 



C(r,L) 



1 



^/2^L 5 / 2 [x(l- x)f/ 2 



(71) 



where x = r/L is the scaling variable. If we identify p 
with m/L, we can identify the expressions for P{p) and 
C(r, L) with the corresponding equilibrium quantities - 
P(n, L) = ±P(p) and G{r, L) = L 2 C(r, L). So, for n > 1 
and r > 1 



dV 2 < F 2 |e- H ( L - r )|Fi > 

J 

<V 1 \e- Hr \V2>e^ +v *\ (68) 



P(n,L) = 



/3 4 L 2 



r 



In 



and 



G(r,L) 



^^[ifl-i)] 3 / 2 



(72) 



(73) 



We see that the scaling forms in these cases are 
similar. A fit to the functional forms shows that these 
equilibrium results reproduce quite well the scaling 
exponents and scaling functions for G(r) and P(n) for 
n > 1 obtained numerically for the noncquilibrium case 
uj = K = 1, as can be seen in Figs. ED and 03 though with 
different values of (3. The correlation function matches 
with /3 ~ 4 while (3 ~ 2.3 describes the probability 
distribution of number data well. However, P(0, L) (and 
thus N occ ) does not agree closely in the two cases. The 
equilibrium case can also be used to shed light on the 
dynamical properties of the nonequilibrium steady state. 
We compared our results for G(t,L) with the density- 
density autocorrelation function in the adiabatic u> — > 
limit. To find the latter, we simulated a surface with 
height field h(r,t) evolving according to KPZ dynamics, 
and evaluated the density using the equilibrium weight 
p(r,t) = e-^^/Z. As shown in [l^, the results with 
(3 = 4 agree with the autocorrelation function in the 
nonequilibrium system, apart from a numerical factor. 

It is surprising that results in this equilibrium limit 
describes the non-equilibrium state so well. In the 
non-equilibrium case, the driving force behind particle 
motion and clustering is the surface fluctuation fluctu- 
ation while the equilibrium case, it is the temperature. 
The common feature in both the cases is the surface 
terrain. Thus, in some region of parameter space the 
surface motion mimics temperature and causes the 
particles to redistribute in a certain way. Why the 
equivalent temperature for various quantities is different 
is not clear and deserves further study. 
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V. FUTURE WORK 

In this paper, we have described our results on the 
problem of particles sliding on a KPZ surface, with 
both the surface and the particles moving in the same 
direction, corresponding to the case of particle advection 
by a noisy Burgers flow. We see that in the steady 
state, the two-point density-density correlation function 
diverges near the origin as a function of distance scaled 
with the system size. This is an indicator of strong 
clustering of particles and the defining characteristic 
of a new kind of state - the strong clustering state (SCS). 

Questions arise about the robustness of the strong clus- 
tering state - Does clustering survive in the case of anti- 
advection where the surface and particles move in oppo- 
site directions to each other? What happens if we change 
the symmetry properties of the driving field, and have 
driving by an Edwards- Wilkinson (EW) surface instead 
of the KPZ surface? Does the phenomenon survive in 
higher dimensions? These questions will be addressed in 
a subsequent paper 0], where it will be shown that the 
steady state is of the SCS kind in all these cases, even 
though the degree of clustering differs from one case to 
another. 
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